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Perturbative two- and three-loop coefficients from large /3 Monte Carlo 
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Perturbative coefficients for Wilson loops and the static quark self-energy are extracted from Monte Carlo 
simulations at large (3 on finite volumes, where all the lattice momenta are large. The Monte Carlo results are in 
excellent agreement with perturbation theory through second order. New results for third order coefficients are 
reported. Twisted boundary conditions are used to eliminate zero modes and to suppress Z3 tunneling. 



1. INTRODUCTION 

Simulations using highly improved lattice ac- 
tions have become commonplace in recent years. 
Higher order perturbative calculations for these 
actions present a major challenge however, which 
must be overcome in order to obtain precision re- 
sults for a number of observables. 

A simple alternative to analytic perturbation 
theory, proposed in Ref. |l|, is to directly mea- 
sure observables in Monte Carlo simulations at 
very large /3 and to fit the coefficients of the per- 
turbative expansion to the results. This produces 
coefficients with both statistical and truncation 
errors, so it is not a complete substitute for con- 
ventional perturbation theory. However, as we 
show here, it produces estimates of high order co- 
efficients with far less effort than a conventional 
calculation. This method was shown to reproduce 
analytical results for the one-loop mass renormal- 
ization for Wilson fcrmions, and the one- loop ad- 
ditive energy for NRQCD fermions Q. In ad- 
dition, preliminary estimates of some third-order 
Wilson loop coefficients have been made ||^. An 
extension of this technique to background field 
calculations has also been considered Q. 

In this work we make a much more extensive 
analysis of this method. We reproduce known 
second-order perturbative coefficients for Wilson 
loops and the static quark self-energy to high ac- 
curacy. New predictions for third-order coeffi- 
cients for these observables are made, and results 



of extensive systematic studies are reported. 

2. WILSON LOOPS 

To begin with, we present results of simula- 
tions of Wilson loops for the Wilson gauge-field 
action on 16^ lattices. Periodic boundary condi- 
tions were used here in order to make a direct 
comparison with the second order perturbative 
coefficients calculated by Heller and Karsch |J] on 
the same lattice volume. Simulations were done 
at nine couplings, from /? « 9 to /? « 50. 

We analyze the logarithm of the RxT Wilson 
loop 
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scaled such that the coefficients for large loops 
approach those of the static self-energy. We use a 
renormalized coupling ap determined from mea- 
sured values of the plaquette ||,^ 

4:7T 

-InW/^i^ = — ap(3.41/a) [1 - 1.185 ap] . 
o 

This eliminates large renormalizations of the bare 
lattice coupling ao = 3/(27r/3), as can be seen by 
relating ap to ao, using the third-order expansion 
of the plaquette given in Ref. m : 

(1) 



ap(3.41/a)wao +4.558 ag 



28.499 ag. 



The renormalized coupling is evaluated at scales 
g^ rp determined by the procedure of Ref. B . It is 



also possible to perform the fitting procedure us- 
ing the bare coupling constant, ao, but this leads 
to much larger higher order coefficients and trun- 
cation errors because of the large coefficients in 
Eqn. |l|. These lead in turn to poorer fits with 
larger x^'s. 

It is very straightforward to fit the leading per- 
turbative coefficients to the Monte Carlo data. 
Monte Carlo results for the 5x5 Wilson loop are 
presented in Fig. m, in terms of the quantity 

Ki = - InWj^jf/ {2{R + r)ap(g|j,T)) , 

after the known first-order contribution due to 
zero modes [B| is subtracted from the data. The 



are shown in Fig. || for the 5x5 loop, with the 
results of a fit to K3 ~ C3 + Ciap. 
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Figure 1. Monte Carlo results for k\ for the 5x5 
Wilson loop, with best fit. 



dashed line in Fig. |l| 



IS the result of a fit to 
Ki « ci -h C2ap -t- c^ap] the filled square shows 
the analytical value of ci , while the curvature in 
the results for ki demonstrates a signal for C3. 

The measured values of ci^2 for various small 
loops are shown in Table |^. They are in excellent 
agreement with perturbation theory. 

We also make new predictions for the third- 
order coefficients C3. To improve the accuracy of 
these results, it is helpful to subtract the analyt- 
ically known first- and second-order perturbative 
contributions from the data. Monte Carlo results 
for the residual 
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Figure 2. Monte Carlo results for K3 for the 5x5 
Wilson loop, with best fit. 



The data are not accurate enough to resolve C4, 
although the best fit errors indicate that C4 is of 
the same order as C3 for these Wilson loops. To 
obtain errors on the C3 that take into account the 
effects of truncation in a^, we include the C4 term 
in our fits. In the results for C3 given in Table 0, 
the quoted errors come from the fit including C4, 
which effectively includes an uncertainty due to 
the possible values of C4 allowed by the data. 

We used periodic boundary conditions in these 
tests on small Wilson loops in order to be able 
to compare directly with Rcf. Q. The known 
leading effects of the resulting zero modes were 
subtracted by hand. We have not corrected for 
the contribution from zero modes beyond first or- 
der; these effects should be very small for small 
Wilson loops. In fact there is no visible effect of 
zero modes beyond first order, within statistical 
errors. (This would show up in singular behavior 
in K3 at small ap^. 

3. STATIC QUARK SELF-ENERGY 

To calculate quantities more complicated than 
Wilson loops, we need to solve some additional 
problems. In general, the effects of zero modes are 



Table 1 

Perturbative coefficients C1.2 of Wilson loops from Monte Carlo simulations and perturbation theory. 



Loop 
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1r,t 


1 X 2 


1.2037(1) 


1.2039 


-1.248(08) 


-1.257 


3.07 


1 X 3 


1.2586(2) 


1.2589 


-1.186(09) 


-1.195 


3.01 


2x2 


1.4336(2) 


1.4338 


-1.314(08) 


-1.320 


2.65 


3x3 


1.6088(3) 


1.6089 


-1.205(12) 


-1.204 


2.46 


4x4 


1.7065(4) 


1.7066 


-1.202(16) 


-1.198 


2.30 


5x5 


1.7692(6) 


1.7690 


-1.191(36) 


-1.166 


2.23 



Table 2 

Monte Carlo results for third order coefficients C3 

of Wilson loops. 



Loop 
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1 X 2 


0.40(08) 


3.07 


1 X 3 


0.63(10) 


3.01 


1 X 4 


0.82(10) 


2.97 


1 X 5 


0.92(10) 


2.95 


2x2 


1.38(09) 


2.65 


2x3 


1.89(11) 


2.56 


2x4 


2.26(11) 


2.49 


2x5 


2.51(11) 


2.46 


3x3 


2.49(12) 


2.46 


3x4 


2.94(11) 


2.38 


3x5 


3.23(13) 


2.35 


4x4 


3.35(14) 


2.30 


4x5 


3.67(19) 


2.27 


5x5 


3.74(33) 


2.23 



unknown and not small, so wc must remove them 
from the first. Further, most quantities are highly 
sensitive to the effects of tunneling between the 
various Z3 vacua of S'[/(3), so we must suppress 
this, too. We test our ability to do this by calcu- 
lating the static quark self-energy Eq. 

Zero modes can be eliminated by using twisted 
boundary conditions [^ , and remaining finite size 
effects can be removed by extrapolating results 
from several lattice volumes. 

We measured the gauge-invariant Polykaov line 



on lattices of various volumes L^ x {T ~ L) 

1 ^ 

PiiL) ^ —Y^^eTr [] U,{x). 

Defining Eq{L) = — lnP4(L)/L, we have Eq = 
Eo{L — > cxd). The two- loop expression for Eq can 
be obtained from perturbative results for P4 p 10 

Eo = 2. 1173 ap(1.68/a)-l. 124 ap-KO(a|,). 

The Polyakov line is an order parameter for 
the Z3 degenerate vacua of the gauge theory, and 
has different values in different Z3 phases. To 
make its vacuum expectation value well defined 
and nonzero in a Monte Carlo simulation, one 
can add an external field to the action to pin the 
simulation into one of the Z^ states and take the 
limit of the external field going to zero as the 
volume goes to infinity. 

In order to minimize these nonperturbative ef- 
fects, we start the simulation with all links ini- 
tialized to U/j, = /; tunneling and domain for- 
mation are suppressed by working at sufficiently 
large lattice volumes and couplings (3. We found 
that with periodic boundary conditions, tunnel- 
ings were quite frequent, even at the surprisingly 
large /3's of 50 ~ 100. On the other hand, we find 
that twisted boundary conditions lead to a dra- 
matic suppression of these effects. Twisting in 
two directions strongly suppresses tunneling and 
twisting in three directions virtually eliminates it, 
even on runs of hundreds of thousands of sweeps 
and relatively low j3. Thus, so far, adopting the 
pinning method has been unnecessary. 

Finite-size effects are consistent with pertur- 
bative expectations. We estimate Eo by fitting 



Eo{L) to the perturbative form of the interaction 
between the static quark and its images in the 
box walls: 

Ea{L) = Eo- (ai + a2ln(L))/L. 

Simulations were done with twisted boundary 
conditions on 9 volumes, L = [3,11]. Measure- 
ments were made at 9 couplings on each volume, 
from /3 w 9 to /3 w 60. Results for Eo{L) at 
/3 = 9.5 are shown in Fig. g|; in all cases, the fits 
yield ai — 0{ap) and 02 




0{ap), as expected. 
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Figure 3. Infinite volume extrapolation of the 
static quark self-energy Eo{L). 



Results for ki 



in Fig. H Best fits 
yield cf^ = 2.117(; 



Eo/ap{q* 

to Kl W 



— 1.68) are shown 
ci + C2ap + c^a'j, 
L17(3) and ci'^ = -1.18(18), in 
good agreement with perturbation theory. We 
also obtain a new prediction for the third order 
term; using the analytical values for ci and C2, 
the best fit yields 
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5.6(1.8). 



(2) 



As for the loops, we include a fourth order term 
in the fit. This guards against the presence of an 
anomalously large C4, and enlarges the error bar. 
A result for the third order term in the expan- 
sion of Eq in the bare lattice coupling has also 
been reported in this conference, using coupled 
Langevin equations mW . Large coefficients in the 



Figure 4. Monte Carlo results for ki for the static 
quark self-energy, with best fit. 



expansion in terms of aa reflect the need to renor- 
malize the coupling; for example, combining our 



result for c^^ with Eqn. (^, wc find 
Eo «2.1173ao + ll-152a^ + 81(2)ajj. 



(3) 



The absolute errors in Eqns. 2 and H are the same. 
The relative error in Eqn. B seems surprisingly 
small because the size of the coefficient itself has 
been increased by analytic terms from Eqn. m. 

4. ALGORITHMIC INVESTIGATIONS 

These calculations have been done using stan- 
dard pseudo-heat bath algorithms to create the 
gauge fields. The calculation of Eq took a little 
over a month running on about 20 PC's. Clearly 
finding the best simulation algorithms is an es- 
sential element of the method. In this section, 
we report some preliminary investigations into al- 
gorithms specifically appropriate to perturbative 
simulations. 

In principle, microcanonical based methods 
III4 have the potential to outperform ordinary 
heat bath methods since they have a critical scal- 
ing dimension of 1 instead of 2. In practice, for 
QCD they don't achieve this. Correlations be- 
tween nearly degenerate modes last a long time, 
and nonperturbative effects cause averages over 
uncorrelated configurations to add up statisti- 
cally, as l/^/n. If we examine the perturbative 



case carefully, however, we see that this need not 
be the case. 

A single harmonic oscillator evolving under a 
microcanonical algorithm approaches its average 
as 1/n: 



(4') 



1 ■c-^ 1 1 - exp(m(S) 

- 2_^ exp{nnS) = - ^ZTTT- (4) 



n 



m=0 



n 1 — cxp(z(5) 



Therefore, free field theory, which is just a 
collection of independent oscillators, does the 
same. However, the usual microcanonical evo- 
lution equations, 



H = T + S, T : 
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\~^ „2 



2^^* 



pf, S = Si(b^) (5) 



system's evolution through phase space deformed 
only slightly. Therefore we expect the covering 
of configuration space to remain dense in weakly 
coupled systems, and the convergence to the av- 
erage to continue to go as 1/n. It is easy to test 
that for scalar field theory, things work out in 
more or less this way. 

In QCD, however, gauge symmetry creates a 
new complication. The random mass method 
may be adapted to gauge theories in two ways 
It is straightforward to implement gauge covari- 
ant random masses that break most but not all 
mode degeneracy: 



T : 



1 
2^ 



TrH? 



S = S{U) 



(9) 



dpi 
dT 



dS{<i>) 



dxJH 

dT 



Pi 



(6) 



do not make a good simulation algorithm for free 
field theory. Many sets of momentum modes, e.g., 
(1, 0, 0, 0) and (0, 1, 0, 0), have same the frequen- 
cies and never decorrelate, so that configuration 
space never gets completely covered. 

Since the "momenta" in microcanonical updat- 
ing of quantum field theories are fictitious, this 
problem and be removed by introducing a "non- 
degenerate microcanonical method" using ran- 
domized fictitious masses to,-: 



H = T + S, 



T : 



dpi 

dT 






lE 



dT 



EL 

rrii 



Pj_ 



5 = ^(00,(7) 



(8) 



ai < rrii < a2, where ai and 02 are of order 1, 
say 0.8 and 1.2. Mathematical classical mechan- 
ics often starts from periodic system system with 
nondegenerate frequencies. Such systems have 
more tractable mathematics, for example no ze- 
ros in perturbation theory denominators. With 
these equations, the frequencies of free field the- 
ory are nondegenerate and configuration space is 
covered densely. The KAM theorem then tells 
us that for small perturbations of nondegener- 
ate unperturbed system, dense tori covered by 



dH, 
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dT 
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(10) 



Mode degeneracy may be broken completely by 
color breaking random masses: 



T = y^^ s = s{u) 

■^—^ TO? 



.dH, 



dSjU) 



dlh 

dT 






(11) 



(12) 



However, they also break gauge invariance. This 
necessitates gauge fixing, which is more work and 
may introduce extra correlations. 

We have investigated the first of these meth- 
ods. New complications are present in gauge 
theories. Free field behavior can only be ap- 
proached as /? — > cx), Up -^ 1. But at very 
high /3, round off becomes important, which de- 
stroys the anticorrelations leading to 1/n behav- 
ior. At small (3, Z3 tunneling occurs. Nonper- 
turbative effects can (and do) destroy the anti- 
correlations, restoring 1/y/n statistical behavior. 
Hence, twisted periodic boundary conditions are 
required. At moderate (3, we find at least par- 
tial success with the plaquette. Fig. || shows the 
error obtained on the plaquette using standard 
binning procedures from the heat bath algorithm 
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Figure 5. Statistical error on the plaquette ob- 
tained with standard binning procedure from the 
heat bath algorithm (canonical) and the nonde- 
generate microcanonical algorithm (microcanon- 
ical). 



(labeled canonical) and the nondegenerate mi- 
crocanonical algorithm (labeled microcanonical). 
When anticorrelations as in Eqn. H are present, 
standard formulae, which assume l/v^ behav- 
ior, should produce errors that drop as l/\Ai 
in bin size, as observed. The results obtained 
are Wi,i = 0.900346(1) with the microcanoni- 
cal method and 0.900271(23) with the canonical 
method for equal ensemble sizes. We have not 
extrapolated to infinite volume, which we believe 
accounts for the small discrepancy. 

The plaquette is speeded up beautifully, achiev- 
ing the full 1/n approach to its average even with 
gauge covariant random masses. However, the 
plaquette turns out to be a special case, presum- 
ably because it is also the action of the theory. 
For other loops, the great speed up is not ob- 
served, which is disappointing but not entirely 
unexpected in view of the remaining correlations. 
Full breaking of mode degeneracy appears neces- 
sary for most quantities. We are encouraged by 



the results with the plaquette, but more work is 
necessary to make this a general method. 

5. SUMMARY 

In summary, results presented here demon- 
strate that second order perturbative coefficients 
are readily accessible in Monte Carlo simulations 
at large (3. Reasonable estimates of third-order 
terms can also be made, especially when analytic 
results for lower order terms are available. We are 
currently working on gauge-fixed quark propaga- 
tors, in order to measure second order mass renor- 
malizations for NRQCD and Fermilab fermions. 

We thank Urs Heller for generously providing 
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This work was supported in part by the National 
Science Foundation, the Department of Energy, 
and the National Science and Engineering Re- 
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